# 稳定性算法
import math

from seaman.core import JsonUtil

from algs.common import RockFormationData


class ShaftResult:
    sigma_xx: float
    sigma_yy: float
    sigma_zz: float
    tau_xy: float
    tau_yz: float
    tau_xz: float


# 井壁围岩应力
def get_shaft_lining_result(data: RockFormationData) -> ShaftResult:
    result = ShaftResult()
    # 求解sigma_xx
    sigma_xx = (data.max_horizontal_ground_stress * math.pow(math.cos(math.radians(data.azimuth)), 2)
                + data.min_horizontal_ground_stress * math.pow(math.sin(math.radians(data.azimuth)), 2)) * math.pow(math.cos(math.radians(data.inclination)), 2) + data.overburden_pressure * math.pow(math.sin(math.radians(data.inclination)), 2)
    result.sigma_xx = sigma_xx
    # 求解sigma_yy
    sigma_yy = data.max_horizontal_ground_stress * math.pow(math.sin(math.radians(data.azimuth)), 2) + data.min_horizontal_ground_stress * math.pow(math.cos(math.radians(data.azimuth)), 2)
    result.sigma_yy = sigma_yy
    # 求解sigma_zz
    sigma_zz = (data.max_horizontal_ground_stress * math.pow(math.cos(math.radians(data.azimuth)), 2)
                + data.min_horizontal_ground_stress * math.pow(math.sin(math.radians(data.azimuth)), 2)) * math.pow(math.sin(math.radians(data.inclination)), 2) + data.overburden_pressure * math.pow(math.cos(math.radians(data.azimuth)), 2)
    result.sigma_zz = sigma_zz
    # 求解tau_xy
    tau_xy = (data.min_horizontal_ground_stress - data.max_horizontal_ground_stress) * 0.5 * math.sin(math.radians(2 * data.azimuth)) * math.cos(math.radians(data.inclination))
    result.tau_xy = tau_xy
    # 求解tau_yz
    tau_yz = (data.min_horizontal_ground_stress - data.max_horizontal_ground_stress) * 0.5 * math.sin(math.radians(2 * data.azimuth)) * math.sin(math.radians(data.inclination))
    result.tau_yz = tau_yz
    # 求解tau_xz
    tau_xz = (data.max_horizontal_ground_stress * math.pow(math.cos(math.radians(data.azimuth)), 2)
              + data.min_horizontal_ground_stress * math.pow(math.sin(math.radians(data.azimuth)), 2)
              - data.overburden_pressure) * 0.5 * math.sin(math.radians(2 * data.inclination))
    result.tau_xz = tau_xz
    return result


class PressureResult:
    pr: float
    m: float

# 孔隙压力
def get_pore_pressure_result(data: RockFormationData) -> PressureResult:
    result = PressureResult()

    # 求解参数 m
    m = 2.12 * math.pow(10, -3) * ((data.fluid_loss * data.drilling_fluid_viscosity)/(data.permeability * data.rock_thickness))
    result.m = m
    # 求解有效孔隙压力 Pr
    # pr = data.in_situ_pore_pressure + m * (((data.conductivity_coefficient * data.drilling_fluid_seepage_time) / math.pow(data.borehole_radius, 2)) - math.log(data.elastic_storage_ratio, 10) + 0.9077)
    pr = data.in_situ_pore_pressure + m * ((data.conductivity_coefficient + math.pow(m, data.drilling_fluid_seepage_time)) / math.pow(data.borehole_radius, 2) - math.log(data.elastic_storage_ratio, 10) + 0.9077)
    result.pr = pr
    return result


class RoundResult:
    z: float
    sigma_r: float
    sigma_theta: float
    sigma_z: float
    tau_theta_z: float


# 井周角求解
def get_round_angle_result(data: RockFormationData) -> RoundResult:
    result = RoundResult()
    temp_one: ShaftResult = get_shaft_lining_result(data)
    sigma_xx = temp_one.sigma_xx
    sigma_yy = temp_one.sigma_yy
    sigma_zz = temp_one.sigma_zz
    tau_xy = temp_one.tau_xy
    tau_yz = temp_one.tau_yz
    tau_xz = temp_one.tau_xz
    temp_two: PressureResult = get_pore_pressure_result(data)
    pr = temp_two.pr
    # pr = 19.751685
    # theta
    theta: float
    # theta = 90
    # pm
    pm: float
    # pm = 1
    # Z
    z = 2 * (sigma_xx - sigma_yy) * math.cos(math.radians(2 * theta)) + 4 * tau_xy * math.sin(math.radians(2 * theta))
    result.z = z

    # sigma_r
    sigma_r = pm - data.osmotic_coefficient * data.void_ratio * (pm - pr)
    result.sigma_r = sigma_r
    # sigma_theta
    sigma_theta = sigma_xx + sigma_yy - z - pm + data.osmotic_coefficient * (pm - pr) * (((data.formation_stress_factor * (1 - 2 * data.poisson_ratio)) / (1 - data.poisson_ratio)) - data.void_ratio)
    result.sigma_theta = sigma_theta
    # sigma_z
    sigma_z = sigma_zz - data.poisson_ratio * z + data.osmotic_coefficient * (pm - pr)*(((data.formation_stress_factor * (1 - 2 * data.poisson_ratio)) / 2 / (1 - data.poisson_ratio)) - data.void_ratio)
    result.sigma_z = sigma_z
    # tau_theta_z
    tau_theta_z = 2 * tau_yz * math.cos(math.radians(theta)) - 2 * tau_xz * math.sin(math.radians(theta))
    result.tau_theta_z = tau_theta_z
    # tau_rz
    tau_rz = 0
    # tau_r_theta
    tau_r_theta = 0

    return result


class EffectiveResult:
    sigma_i: float
    sigma_j: float
    sigma_k: float

# 井周有效应力
def get_well_effective_stress_result(data: RockFormationData) -> EffectiveResult:
    result = EffectiveResult()
    temp_one = get_pore_pressure_result(data)
    pr = temp_one.pr
    # pr = 19.751685
    temp_two = get_round_angle_result(data)
    sigma_r = temp_two.sigma_r
    sigma_theta = temp_two.sigma_theta
    # sigma_theta = 67.03039908
    sigma_z = temp_two.sigma_z
    tau_theta_z = temp_two.tau_theta_z

    # sigma_i
    sigma_i = sigma_r - data.formation_stress_factor * pr
    result.sigma_i = sigma_i
    # 定义右部分子计算元
    factor_a = math.pow(sigma_theta - sigma_z, 2) + 4 * math.pow(tau_theta_z, 2)
    # sigma_j
    sigma_j = (sigma_theta + sigma_z) * 0.5 + math.pow(factor_a, 0.5) * 0.5 - data.formation_stress_factor * pr
    result.sigma_j = sigma_j
    # sigma_k
    sigma_k = (sigma_theta + sigma_z) * 0.5 - math.pow(factor_a, 0.5) * 0.5 - data.formation_stress_factor * pr
    result.sigma_k = sigma_k

    return result


class DpResult:
    alpha: float
    k: float
    j2: float
    i1: float
    result_y: float

# 失稳准则： D.P准则
def get_dp_result(data: RockFormationData) -> DpResult:
    result = DpResult()
    temp = get_well_effective_stress_result(data)
    sigma_i = temp.sigma_i
    sigma_j = temp.sigma_j
    sigma_k = temp.sigma_k

    # 根号下的计算元 factor_a
    factor_a = math.sin(math.radians(data.internal_friction_angle))
    # alpha
    alpha = (math.pow(3, 0.5) * factor_a) / (3 * math.pow((3 + math.pow(factor_a, 2)), 0.5))
    result.alpha = alpha
    # K
    k = (math.pow(3, 0.5) * data.cohesion * math.cos(math.radians(data.internal_friction_angle))) / math.pow((3 + math.pow(factor_a, 2)), 0.5)
    result.k = k
    # j2
    j2 = (math.pow(sigma_i - sigma_j, 2) + math.pow(sigma_j - sigma_k, 2) + math.pow(sigma_k - sigma_i, 2)) / 6
    result.j2 = j2
    # i1
    i1 = sigma_i + sigma_j + sigma_k
    result.i1 = i1
    # result_y
    result_y = math.pow(j2, 0.5) - alpha * i1 - k
    result.result_y = result_y
    # if result_y > 0:
    #     return "井壁失稳"
    # elif result_y < 0:
    #     return "井壁稳定"
    return result


if __name__ == '__main__':
    data = RockFormationData()
    data.depth = 2000
    data.inclination = 30
    data.azimuth = 0
    data.cohesion = 2
    data.internal_friction_angle = 34
    data.formation_stress_factor = 0.7
    data.poisson_ratio = 0.2
    data.max_horizontal_ground_stress = 36.6
    data.min_horizontal_ground_stress = 35.6
    data.overburden_pressure = 41
    data.in_situ_pore_pressure = 19.6
    data.borehole_radius = 10.8
    data.void_ratio = 0.02
    data.fluid_loss = 1
    data.drilling_fluid_viscosity = 50
    data.elastic_storage_ratio = 0.8
    data.conductivity_coefficient = 80000
    data.drilling_fluid_seepage_time = 3600
    data.permeability = 120
    data.rock_thickness = 4
    data.osmotic_coefficient = 1

    shaft_result = get_shaft_lining_result(data)
    pressure_result = get_pore_pressure_result(data)
    round_result = get_round_angle_result(data)
    well_result = get_well_effective_stress_result(data)
    dp_result = get_dp_result(data)
    print('===当前压力值为===')
    print(JsonUtil.obj2str(dp_result))